%     AERO LOAD - программа для нагружения агрегатов самолета распределенной нагрузкой
%     Copyright (C) 2018  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

more off;
clc;

% #######################################################################################################
% ##                                                                                                   ##
% ##                                        НАСТРОЙКИ ПРОГРАММЫ                                        ##
% ##                                                                                                   ##
% #######################################################################################################

% Для включения возможности вывода информации на русском языке при выполнении программы в среде GNU Octave
% необходимо добавить в файл C:\Octave\Octave-4.4.0\share\octave\site\m\startup\octaverc  следующую строку:
% system("chcp 1251",1);

% --------------------------------------- Параметры ввода-вывода ----------------------------------------

Work_Directory                 = 'M:\Patran_work'; % Рабочая папка
Need_Reading_BDF_File          = 'TRUE'; % TRUE | FALSE     % TRUE - выполнить новое чтение bdf файла; FALSE - использовать базу, содержащуюся в памяти программы
Input_BDF_name                 = 'example'; % Имя bdf файла

% -------------------------------------- ЗАДАНИЕ БАЗОВОЙ ПЛОСКОСТИ --------------------------------------

% Базовая плоскость - это та плоскость, перпендикулярно которой должна действовать нагрузка

% Способ задания базовой плоскости
Basic_Plane_Definition = 'THREE_POINTS'; % POINT_NORMAL | THREE_POINTS

% Если Basic_Plane_Definition = POINT_NORMAL (Плоскость задается по точке Z0 и нормали Z0 -> Z1)
Z0 = [207; -32; 161];
Z1 = [204.676; -111.973; 204.359];

% Если Basic_Plane_Definition = THREE_POINTS (Плоскость задается по трем точкам D1, D2, D3; нормаль определяется последовательностью обхода D1 -> D2 -> D3)
D1 = [0; 0; 0];
D2 = [0; -400; 0];
D3 = [2000; 0; 0];

% -------------------------------- ЗАДАНИЕ УГЛОВЫХ ТОЧЕК ПОВЕРХНОСТИ ХОРД -------------------------------

% Угловые точки поверхности хорд будут спроецированы на базовую плоскость
% На проекции поверхности хорд на базовую плоскость будет воссоздано поле распределенной нагрузки

% A1 -> A2 - направление размаха по передней кромке
A1 = [0; 0; 100];
A2 = [2000; 0; 100];

% A3 -> A4 - направление размаха по задней кромке
A3 = [0; -400; 100];
A4 = [2000; -200; 100];

% ---------------------------------------- ЗАДАНИЕ СУММАРНОЙ СИЛЫ ---------------------------------------

% Суммарная сила, действующая по нормали к базовой плоскости
Pn = 3000;

% ---------------------- ЗАДАНИЕ ЗАКОНА РАСПРЕДЕЛЕНИЯ ПОГОННОЙ НАГРУЗКИ ПО РАЗМАХУ ----------------------

% Закон распределения погонной нагрузки по размаху
Span_Load_Law = 'TRAPEZOID'; % CHORD_PROPORTIONAL | TRAPEZOID

    % CHORD_PROPORTIONAL - распределение по трапеции, построенной пропорционально длинам корневой и концевой хорд. Длины корневой и концевой хорд - это длины спроецированных на базовую плоскость отрезков A1A3 и A2A4
    % TRAPEZOID - распределение по трапеции с заданным центром давления

    % Если Span_Load_Law = TRAPEZOID
    Span_PC_Type = 'RATIO'; % RATIO | LOCATION (способ задания центра давления)

        % Если Span_PC_Type = RATIO (Центр давления задается как параметр, 0.333(3) <= Span_PC_Ratio <= 0.666(6) )
        Span_PC_Ratio = 0.45;

        % Если Span_PC_Type = LOCATION (Центр давления задается координатами точки. Эта точка проецируется на базовую плоскость, после чего за центр давления принимается ее параметрическая координата по размаху)
        Span_PC_Location = [1346.729; 248.456; 739.358];

% --------------------------- ЗАДАНИЕ ЗАКОНА РАСПРЕДЕЛЕНИЯ ДАВЛЕНИЯ ПО ХОРДАМ ---------------------------

% Закон распределения давления по хордам
Chord_Load_Law = 'TRAPEZOID'; % TRIANGLE | TRAPEZOID

    % TRIANGLE - распределение по треугольнику (центр давления 0.333(3) )
    % TRAPEZOID - распределение по трапеции с заданным центром давления

    % Если Chord_Load_Law = TRAPEZOID
    Chord_PC_Type = 'RATIO'; % RATIO | LOCATION (способ задания центра давления)

        % Если Chord_PC_Type = RATIO (Центр давления задается как параметр, 0.333(3) <= Chord_PC_Ratio <= 0.666(6) )
        Chord_PC_Ratio = 0.4;

        % Если Chord_PC_Type = LOCATION (Центр давления задается координатами точки. Эта точка проецируется на базовую плоскость, после чего за центр давления принимается ее параметрическая координата по хорде)
        Chord_PC_Location = [1346.729; 248.456; 739.358];

% ----------------------------------- НАСТРОЙКИ СУММИРОВАНИЯ НАГРУЗОК -----------------------------------

% Суммирование нагрузок в программе Aero Load используется в целях проверки фактической нагрузки, приходящей на элементы

Summation_Point = [-100; -50; 0]; % Точка суммирования (влияет на результирующий момент)

% Прямоугольная система координат для суммирования нагрузок:
Coord_Origin = [0; 0; 0];    % Начало системы координат
Point_On_Axis_3 = [0; 0; 1];    % Точка на оси Z
Point_On_Plane_13 = [1; 0; 0];    % Точка в плоскости XZ

% ------------------------------------ НАСТРОЙКИ КОРРЕКЦИИ НАГРУЗОК -------------------------------------

Correction_Coefficient = 1; % Поправочный коэффициент, на который будет умножена приходящая на элементы нагрузка

% ---------------------------------------- ПАРАМЕТРЫ СОХРАНЕНИЯ -----------------------------------------

Vector_Field_Component = 3; % 1 | 2 | 3    % Номер компоненты векторного поля для CID Distributed Load, в который будет записываться нагрузка. Остальные компоненты векторного поля будут равны нулю
Output_File_Format = 'CSV'; % SES | CSV | BOTH | OFF    % Формат выходного файла (SES, CSV, или оба сразу), либо отключение (OFF) сохранения поля нагрузок в файл (полезно при выполнении коррекции нагрузки)

    % Если Output_File_Format = SES
    Field_Name = 'Field_P_case_A'; % Имя поля, которое будет создано в Patran

% -------------------------------------------------------------------------------------------------------

% #######################################################################################################
% ##                                                                                                   ##
% ##                                     КОНЕЦ НАСТРОЕК ПРОГРАММЫ                                      ##
% ##                                                                                                   ##
% #######################################################################################################

% Очистка памяти
if (strcmp(Need_Reading_BDF_File, 'TRUE'))

    clear ErrCode;
    
    clear DB_NODES;
    clear DB_ELEMENTS;
    clear DB_MPC;
    clear DB_MATERIALS;
    clear DB_PSETS;
    clear DB_COORD;
    clear DB_FIELDS;
    
    clear AT_NODES;
    clear AT_ELEMENTS;
    clear AT_MPC;
    clear AT_MATERIALS;
    clear AT_PSETS;
    clear AT_COORD;
    clear AT_FIELDS;
    
    clear db_Next_Field_Number;
    clear BDF_Entity_Read_Settings;

end

% -----------------------------------------------

global software_version; %                      1
global software_date; %                         1
global software_year; %                         1

global ErrCode; %                               1
global Progress_Bar_Segments; %                 1
global Progress_Bar_Delay; %                    1
global Progress_Bar_Elm_Threshold; %            1

global MAT2_orthotropic_tolerance_multiplier; % 1
global MAT9_orthotropic_tolerance_multiplier; % 1
global MAT2_C_zero_det_tolerance; %             1
global MAT9_C_zero_det_tolerance; %             1
global Laminate_offset_round_tolerance; %       1
global Orientation_Vector_zero_tolerance; %     1
global Global_equivalence_tolerance; %          1
global Field_Zero_Record_Tolerance; %           1
global Field_Nonzero_Record_Tolerance; %        1

global DB_NODES; %                              1
global DB_ELEMENTS; %                           1
global DB_MPC; %                                1
global DB_MATERIALS; %                          1
global DB_PSETS; %                              1
global DB_COORD; %                              1
global DB_FIELDS; %                             1

global AT_NODES; %                              1
global AT_ELEMENTS; %                           1
global AT_MPC; %                                1
global AT_MATERIALS; %                          1
global AT_PSETS; %                              1
global AT_COORD; %                              1
global AT_FIELDS; %                             1

global AT_Page_Prealloc_Size; %                 1
global db_Next_Field_Number; %                  1
global BDF_Entity_Read_Settings; %              1

global CID_Vector_zero_tolerance;
global db_Next_Coord_Number;

% -----------------------------------------------

ErrCode = 0;

if (strcmp(Need_Reading_BDF_File, 'TRUE'))
    
    software_version = [];
    software_date = [];
    software_year = [];

    Progress_Bar_Segments = [];
    Progress_Bar_Delay = [];
    Progress_Bar_Elm_Threshold = [];

    MAT2_orthotropic_tolerance_multiplier = [];
    MAT9_orthotropic_tolerance_multiplier = [];
    MAT2_C_zero_det_tolerance = [];
    MAT9_C_zero_det_tolerance = [];
    Laminate_offset_round_tolerance = [];
    Orientation_Vector_zero_tolerance = [];
    Global_equivalence_tolerance = [];
    Field_Zero_Record_Tolerance = [];
    Field_Nonzero_Record_Tolerance = [];
    
    CID_Vector_zero_tolerance = [];

    DB_FIELDS = [];
    AT_FIELDS = [];

    AT_Page_Prealloc_Size = [];
    db_Next_Field_Number = 1;
    db_Next_Coord_Number = [];
    BDF_Entity_Read_Settings = [];

end

% -----------------------------------------------

% Отсчет времени выполнения программы
tic;

% Инициализация базовых переменных
Init_Variables ();

% Инициализация допусков
Init_Tolerances ();

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Печать приветствия
Print_Greeting ();

% Корректировка Work_Directory

if (Work_Directory(length(Work_Directory)) == '\')
    Work_Directory = Work_Directory(1:length(Work_Directory)-1);
end

% Корректировка имени bdf файла

if ((length(Input_BDF_name) > 4) && strcmp(Input_BDF_name(length(Input_BDF_name)-3:length(Input_BDF_name)), '.bdf'))
    Input_BDF_name = Input_BDF_name(1:length(Input_BDF_name)-4);
end

% Полный путь к bdf файлу

BDF_File_Path = [Work_Directory '\' Input_BDF_name '.bdf'];

% ======== Этап создания модели ========

if (strcmp(Need_Reading_BDF_File, 'TRUE'))

    fprintf (1, 'СОЗДАНИЕ МОДЕЛИ\n\n');

    % Чтение bdf файла

    [DB_NODES, ELEMENTS, MPC, MATERIALS, PSETS, COORD] = Read_BDF(BDF_File_Path);

    if (ErrCode < 0)
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
    end

    % Создание адресных таблиц

    [AT_COORD, AT_NODES, AT_MPC, AT_MATERIALS, AT_ELEMENTS, AT_PSETS] = Create_Address_Tables(COORD, MPC, MATERIALS, ELEMENTS, PSETS);

     if (ErrCode < 0)
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
    end

    % --- Проверка объектов ---

    % Проверка узлов

    if (~isempty(DB_NODES))

        Check_Nodes();

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end

    end

    % Проверка MPC

    if (~isempty(MPC))

        Check_MPC(MPC);

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end

    end

    % Проверка материалов

    if (~isempty(MATERIALS))

        Check_Materials(MATERIALS);

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end

    end

    % Проверка property sets

    if (~isempty(PSETS))

        Check_Psets(PSETS, MATERIALS, COORD);

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end

    end

    % Проверка элементов

    if (~isempty(ELEMENTS))

        Check_Elements(ELEMENTS, PSETS, COORD);

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end

    end

    % --- Обработка объектов и создание баз данных ---

    % Обработка систем координат

    if (~isempty(COORD))

        DB_COORD = Create_DB_COORD(COORD);

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end

        clear COORD;

    end

    % Обработка узлов

    if (~isempty(DB_NODES))

        Fill_DB_NODES_Locations();

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end

        clear NODES;

    end

    % Обработка MPC

    if (~isempty(MPC))

        DB_MPC = Create_DB_MPC(MPC);

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end

        clear MPC;

    end

    % Обработка материалов

    if (~isempty(MATERIALS))

        DB_MATERIALS = Create_DB_MATERIALS(MATERIALS);

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end

        clear MATERIALS;

    end

    % Обработка property sets

    if (~isempty(PSETS))

        DB_PSETS = Create_DB_PSETS(PSETS);

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end

        clear PSETS;

    end

    % Обработка элементов

    if (~isempty(ELEMENTS))

        DB_ELEMENTS = Create_DB_ELEMENTS(ELEMENTS);

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end

        clear ELEMENTS;

    end

    % --- Вывод статистики ---

    fprintf (1, '\nСтатистика модели:\n\n');

    if (strcmp(BDF_Entity_Read_Settings.Read_NODES, 'ON'))
        fprintf (1, '               Узлов: %i\n', size(DB_NODES,1));
    end

    if (strcmp(BDF_Entity_Read_Settings.Read_PSETS_ELEMENTS, 'ON'))
        fprintf (1, '           Элементов: %i\n', size(DB_ELEMENTS,1));
    end

    if (strcmp(BDF_Entity_Read_Settings.Read_MPC, 'ON'))
        fprintf (1, '                 MPC: %i\n', size(DB_MPC,1));
    end

    if (strcmp(BDF_Entity_Read_Settings.Read_MATERIALS, 'ON'))
        fprintf (1, '          Материалов: %i\n', size(DB_MATERIALS,1));
    end

    if (strcmp(BDF_Entity_Read_Settings.Read_PSETS_ELEMENTS, 'ON'))
        fprintf (1, '       Property sets: %i\n', size(DB_PSETS,1));
    end

    if (strcmp(BDF_Entity_Read_Settings.Read_COORD, 'ON'))
        fprintf (1, '    Систем координат: %i\n', size(DB_COORD,1));
    end

    Print_Delimiter();

end


% ------------------------------------------------ ФОРМАТ ГЛОБАЛЬНЫХ БАЗ ДАННЫХ ------------------------------------------------

% DB_NODES | .ID | .Ref_Location(3,1) | .Location(3,1) | .Reference_CID | .Analysis_CID |

% DB_ELEMENTS (Tet | Wedge | Hex) | .type | .ID | .Prop_ID | .N_nodes | .Nodes(N_nodes,1) |
% DB_ELEMENTS (Quad | Tria)       | .type | .ID | .Prop_ID | .N_nodes | .Nodes(N_nodes,1) |
% DB_ELEMENTS (Bar)               | .type | .ID | .Prop_ID | .N_nodes | .Nodes(N_nodes,1) |

% DB_MPC (RBAR) | .type | .ID | .ind_grid | .ind_DOF (в СК анализа узла ind) | .dep_grid | .dep_DOF (в СК анализа узла dep) |
% DB_MPC (RBE2) | .type | .ID | .ind_grid | .dep_grid(:,1) | .dep_DOF (в СК анализа узлов dep) |
% DB_MPC (RBE3) | .type | .ID | .ind_data(:,3) | .dep_grid | .dep_DOF (в СК анализа узла dep) |
%            .ind_data(i,:) | weight_i | ind_grid_i | dof_i (в СК анализа узла ind_grid_i) |

% DB_MATERIALS (isotropic)      | .type | .ID | .Patran_name | .E | .nu |
% DB_MATERIALS (2d_orthotropic) | .type | .ID | .Patran_name | .E1 | .E2 | .nu12 | .G12 |
% DB_MATERIALS (2d_anisotropic) | .type | .ID | .Patran_name | .C(3,3) - матрица Гука |                                         ### Сдвиговые деформации инженерные
% DB_MATERIALS (3d_orthotropic) | .type | .ID | .Patran_name | .E1 | .E2 | .E3 | .G12 | .G23 | .G13 | .nu12 | .nu23 | .nu13 |   ### nu12*E2 = nu21*E1 # nu23*E3 = nu32*E2 # nu13*E3 = nu31*E1
% DB_MATERIALS (3d_anisotropic) | .type | .ID | .Patran_name | .C(6,6) - матрица Гука |                                         ### Сдвиговые деформации инженерные # Расшифровка номеров сдвиговых компонент: 4 - 12, 5 - 23, 6 - 13 # nu12*E2 = nu21*E1 # nu23*E3 = nu32*E2 # nu13*E3 = nu31*E1
% DB_MATERIALS (pshell_MAT2_for_transverse_shear) | .type | .ID | .Patran_name | .C(2,2) - матрица Гука |

% DB_PSETS ( PSOLID) | .type | .ID | .Patran_name | .Mat_ID | .Ori_CID (-1, если element coordinate system) |
% DB_PSETS ( PCOMP)  | .type | .ID | .Patran_name | .Patran_laminate_name | .N_plies | .Laminate(:,3) | .Laminate_offset |
%                    | .Ori_Field_ID | .Offset_Field_ID |         DB_FIELDS{Ori_Field_ID}.Value_data(:,2) [Value_data(j,1) - ori_type: (0 - CID | 1 - Angle)] [Value_data(j,2) - ori_data]
%             .Laminate(i,3) | Mat_ID_i | t_i | Theta_i |
% DB_PSETS ( PCOMPG) | .type | .ID | .Patran_name | .Patran_laminate_name | .N_plies | .Laminate(:,4) | .Laminate_offset |
%                    | .Ori_Field_ID | .Offset_Field_ID |         DB_FIELDS{Ori_Field_ID}.Value_data(:,2) [Value_data(j,1) - ori_type: (0 - CID | 1 - Angle)] [Value_data(j,2) - ori_data]
%             .Laminate(i,4) | Mat_ID_i | t_i | Theta_i | Global_ply_ID_i |
% DB_PSETS ( PSHELL) | .type | .ID | .Patran_name | .Mat_ID_membrane | .Mat_ID_bending | .Mat_ID_transverse_shear | .t |
%                    | .Ori_Field_ID | .Offset_Field_ID |         DB_FIELDS{Ori_Field_ID}.Value_data(:,2) [Value_data(j,1) - ori_type: (0 - CID | 1 - Angle)] [Value_data(j,2) - ori_data]
% DB_PSETS (  PBUSH) | .type | .ID | .Patran_name | .K(6,1) | .Ori_CID_Field_ID |
% DB_PSETS (   PROD) | .type | .ID | .Patran_name | .Mat_ID | .Area | .J |
% DB_PSETS (PBMSECT) | .type | .ID | .Patran_name | .Mat_ID | .Cross_section_form | .Section_data |
%                    | .Ori_vector_Field_ID | .Offset_vector_A_Field_ID | .Offset_vector_B_Field_ID | .Pinned_DOFs_A_Field_ID | .Pinned_DOFs_B_Field_ID |     (вектора ориентации и офсетов в глобальной СК)
% DB_PSETS (PBRSECT) | .type | .ID | .Patran_name | .Mat_ID | .Cross_section_form | .Section_data |
%                    | .Ori_vector_Field_ID | .Offset_vector_A_Field_ID | .Offset_vector_B_Field_ID | .Pinned_DOFs_A_Field_ID | .Pinned_DOFs_B_Field_ID |     (вектора ориентации и офсетов в глобальной СК)
% DB_PSETS (   PBAR) | .type | .ID | .Patran_name | .Mat_ID | .Area | .I1 | .I2 | .I12 | .J | .K1 | .K2 |
%                    | .Ori_vector_Field_ID | .Offset_vector_A_Field_ID | .Offset_vector_B_Field_ID | .Pinned_DOFs_A_Field_ID | .Pinned_DOFs_B_Field_ID |     (вектора ориентации и офсетов в глобальной СК)
% DB_PSETS (  PBARL) | .type | .ID | .Patran_name | .Mat_ID | .Shape | .Dim(:,1) |
%                    | .Ori_vector_Field_ID | .Offset_vector_A_Field_ID | .Offset_vector_B_Field_ID | .Pinned_DOFs_A_Field_ID | .Pinned_DOFs_B_Field_ID |     (вектора ориентации и офсетов в глобальной СК)
%             .Shape = (ROD, TUBE, TUBE2, L, I, CHAN, T, BOX, BAR, CROSS, H, T1, I1, CHAN1, Z, CHAN2, T2, BOX1, HEXA, HAT, HAT1, DBOX )
% DB_PSETS (  PBEAM) | .type | .ID | .Patran_name | .Mat_ID | .Area | .I1 | .I2 | .I12 | .J | .K1 | .K2 | .S1 | .S2 | .N1_A | .N2_A | .N1_B | .N2_B |
%                    | .N_intermediate_states | .Intermediate_states_data{N_intermediate_states,1} |
%                    | .Ori_vector_Field_ID | .Offset_vector_A_Field_ID | .Offset_vector_B_Field_ID | .Pinned_DOFs_A_Field_ID | .Pinned_DOFs_B_Field_ID |     (вектора ориентации и офсетов в глобальной СК)
%             .Intermediate_states_data{i} : { .Location | .A | .I1 | .I2 | .I12 | .J }
% DB_PSETS ( PBEAML) | .type | .ID | .Patran_name | .Mat_ID | .Shape | .Dim(:,1) | .N_intermediate_states | .Intermediate_states_data{N_intermediate_states,1} |
%                    | .Ori_vector_Field_ID | .Offset_vector_A_Field_ID | .Offset_vector_B_Field_ID | .Pinned_DOFs_A_Field_ID | .Pinned_DOFs_B_Field_ID |     (вектора ориентации и офсетов в глобальной СК)
%             .Intermediate_states_data{i} : { .Location | .Dim(:,1) }
%             .Shape = (ROD, TUBE, TUBE2, L, I, CHAN, T, BOX, BAR, CROSS, H, T1, I1, CHAN1, Z, CHAN2, T2, BOX1, HEXA, HAT, HAT1, DBOX )
% DB_PSETS ( PBCOMP) | .type | .ID | .Patran_name | .Mat_ID | .Area | .I1 | .I2 | .I12 | .J | .K1 | .K2 | .N1 | .N2 | .SYMOPT | .N_lumped_areas | .Lumped_area_data{N_lumped_areas,1} |
%                    | .Ori_vector_Field_ID | .Offset_vector_A_Field_ID | .Offset_vector_B_Field_ID | .Pinned_DOFs_A_Field_ID | .Pinned_DOFs_B_Field_ID |     (вектора ориентации и офсетов в глобальной СК)
%             .Lumped_area_data{i} : { .Y | .Z | .C | .MID }

% DB_COORD | .ID | .type (Rectangular | Cylindrical) | .e_x (.e_r) | .e_y(.e_t) | .e_z | .origin |

% DB_FIELDS | .ID | .Dim | .Default | .AT{:,1} | .ID_data{:,1} | .Value_data{:,1} |

% ------------------------------------------------------------------------------------------------------------------------------

% -------------------------------------------------- ПЕРЕЧЕНЬ АДРЕСНЫХ ТАБЛИЦ --------------------------------------------------

% AT_NODES(:,1)
% AT_ELEMENTS(:,1)
% AT_MPC(:,1)
% AT_MATERIALS(:,1)
% AT_PSETS(:,1)
% AT_COORD(:,1)
% AT_FIELDS(:,1)

% ------------------------------------------------------------------------------------------------------------------------------

% ----------------------------------------- ФУНКЦИИ ДЛЯ ДОСТУПА К СВОЙСТВАМ ЭЛЕМЕНТОВ ------------------------------------------

%               [Ori_type, Ori_data] = Get_Shell_Orientation(Element_ID)
%                             Offset = Get_Shell_Offset(Element_ID)    (см. описание функции)
%                             Offset = Get_Plate_Offset(Element_ID)    (см. описание функции)
%                            Ori_CID = Get_Bush_Ori_CID(Element_ID)
%                         Ori_vector = Get_Bar_Ori_Vector(Element_ID)
% [Offset_vector_A, Offset_vector_B] = Get_Bar_Offset_Vectors(Element_ID)
%     [Pinned_DOFs_A, Pinned_DOFs_B] = Get_Bar_Pinned_DOFs(Element_ID)

% ------------------------------------------------------------------------------------------------------------------------------


% ##############################################################################################################################
% ##                                                                                                                          ##
% ##   С ЭТОГО МОМЕНТА НА ОСНОВЕ ДАННОГО КОДА МОЖНО СОЗДАВАТЬ ДРУГИЕ ПРОГРАММЫ, ИСПОЛЬЗУЮЩИЕ ДАННЫЕ ИЗ СФОРМИРОВАННОЙ БАЗЫ    ##
% ##                                                                                                                          ##
% ##############################################################################################################################


% ======== Этап нагружения модели ========

fprintf (1, 'НАГРУЖЕНИЕ МОДЕЛИ\n\n');

% Исправление переменных

if (~iscolumn(Z0))
    Z0 = Z0';
end

if (~iscolumn(Z1))
    Z1 = Z1';
end

if (~iscolumn(D1))
    D1 = D1';
end

if (~iscolumn(D2))
    D2 = D2';
end

if (~iscolumn(D3))
    D3 = D3';
end

if (~iscolumn(A1))
    A1 = A1';
end

if (~iscolumn(A2))
    A2 = A2';
end

if (~iscolumn(A3))
    A3 = A3';
end

if (~iscolumn(A4))
    A4 = A4';
end

if (~iscolumn(Span_PC_Location))
    Span_PC_Location = Span_PC_Location';
end

if (~iscolumn(Chord_PC_Location))
    Chord_PC_Location = Chord_PC_Location';
end

if (~iscolumn(Summation_Point))
    Summation_Point = Summation_Point';
end

if (~iscolumn(Coord_Origin))
    Coord_Origin = Coord_Origin';
end

if (~iscolumn(Point_On_Axis_3))
    Point_On_Axis_3 = Point_On_Axis_3';
end

if (~iscolumn(Point_On_Plane_13))
    Point_On_Plane_13 = Point_On_Plane_13';
end

% Проверка памяти

if (strcmp(Need_Reading_BDF_File, 'FALSE'))
    if (isempty(DB_ELEMENTS))
        ErrCode = -1;
        fprintf (1, 'В памяти программы отсутствует загруженная модель! Измените параметр Need_Reading_BDF_File.\n');
    end
end

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Проверка параметра Vector_Field_Component

if (~((Vector_Field_Component == 1) || (Vector_Field_Component == 2) || (Vector_Field_Component == 3)))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Vector_Field_Component [%i] !\n', Vector_Field_Component);
end

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Проверка элементов

fprintf (1, 'Проверка элементов...                           ');

for i=1:size(DB_ELEMENTS,1)
    switch DB_ELEMENTS{i}.type
        case 'Quad'
            if (DB_ELEMENTS{i}.N_nodes ~= 4)
                ErrCode = -1;
                fprintf (1, 'Недопустимое количество узлов [%i] у элемента с номером %i типа %s !\n', DB_ELEMENTS{i}.N_nodes, DB_ELEMENTS{i}.ID, DB_ELEMENTS{i}.type);
                fprintf(1, 'В программе поддерживаются только Quad элементы с 4 узлами.\n');
                break;
            end
        case 'Tria'
            if (DB_ELEMENTS{i}.N_nodes ~= 3)
                ErrCode = -1;
                fprintf (1, 'Недопустимое количество узлов [%i] у элемента с номером %i типа %s !\n', DB_ELEMENTS{i}.N_nodes, DB_ELEMENTS{i}.ID, DB_ELEMENTS{i}.type);
                fprintf(1, 'В программе поддерживаются только Tria элементы с 3 узлами.\n');
                break;
            end
        otherwise
            ErrCode = -1;
            fprintf (1, 'Недопустимый тип элемента с номером %i [%s] !\n', DB_ELEMENTS{i}.ID, DB_ELEMENTS{i}.type);
            fprintf(1, 'В bdf файле должны содержаться только элементы типов Quad и Tria.\n');
            break;
    end
end

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

fprintf (1, 'OK!\n');

% Уравнение базовой плоскости

switch Basic_Plane_Definition
    case 'POINT_NORMAL'
        Basic_Plane = Plane_Equation_Point_Normal(Z0, Z1);
    case 'THREE_POINTS'
        Basic_Plane = Plane_Equation_3_Points(D1, D2, D3);
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Basic_Plane_Definition [%s] !\n', Basic_Plane_Definition);
end

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Вектор действия нагрузки

Load_Vector = [Basic_Plane(1); Basic_Plane(2); Basic_Plane(3)];
Load_Vector = Load_Vector / norm(Load_Vector);

% Проецирование угловых точек поверхности хорд на базовую плоскость

A1 = Project_Point_To_Plane(A1, Basic_Plane);
A2 = Project_Point_To_Plane(A2, Basic_Plane);
A3 = Project_Point_To_Plane(A3, Basic_Plane);
A4 = Project_Point_To_Plane(A4, Basic_Plane);

% Длины корневой и концевой хорд

Root_chord_length = norm(A1 - A3);
End_chord_length = norm(A2 - A4);

% Центр давления по размаху

switch Span_Load_Law
    case 'CHORD_PROPORTIONAL'
        Span_Pressure_Center = 1/3 * (Root_chord_length + 2*End_chord_length) / (Root_chord_length + End_chord_length);
    case 'TRAPEZOID'
        switch Span_PC_Type
            case 'RATIO'
                Span_Pressure_Center = Span_PC_Ratio;
            case 'LOCATION'
                Span_Pressure_Center = Get_Param_Coord(A1, A2, A3, A4, Project_Point_To_Plane(Span_PC_Location, Basic_Plane));
            otherwise
                ErrCode = -1;
                fprintf (1, 'Недопустимое значение параметра Span_PC_Type [%s] !\n', Span_PC_Type);
        end
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Span_Load_Law [%s] !\n', Span_Load_Law);
end

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Центр давления по хорде

switch Chord_Load_Law
    case 'TRIANGLE'
        Chord_Pressure_Center = 1/3;
    case 'TRAPEZOID'
        switch Chord_PC_Type
            case 'RATIO'
                Chord_Pressure_Center = Chord_PC_Ratio;
            case 'LOCATION'
                Chord_Pressure_Center = Get_Param_Coord(A1, A3, A2, A4, Project_Point_To_Plane(Chord_PC_Location, Basic_Plane));
            otherwise
                ErrCode = -1;
                fprintf (1, 'Недопустимое значение параметра Chord_PC_Type [%s] !\n', Chord_PC_Type);
        end
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Chord_Load_Law [%s] !\n', Chord_Load_Law);
end

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Длина размаха, используемая для расчета погонной нагрузки

Span_length = Get_Point_Line_Distance(A2 - (A2 - A4)*Chord_Pressure_Center, A1, A3);

% Проецирование узлов на базовую плоскость

fprintf (1, 'Проецирование элементов на базовую плоскость... ');

PROJECTED_NODES = zeros(size(DB_NODES,1),3);

for i=1:size(DB_NODES,1)
    PROJECTED_NODES(i,:) = Project_Point_To_Plane(DB_NODES{i}.Location, Basic_Plane)';
end

% Выделение памяти для хранения поля нагрузок

LOAD_FIELD = zeros(size(DB_ELEMENTS,1),4);

fprintf (1, 'OK!\n');

% Вычисление нагрузок на элементы

fprintf (1, 'Вычисление нагрузок на элементы...              ');

for i=1:size(DB_ELEMENTS,1)
    
    LOAD_FIELD(i,1) = DB_ELEMENTS{i}.ID;
    
    switch DB_ELEMENTS{i}.N_nodes
        
        case 3
            
            Node_1 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(1))}.Location;
            Node_2 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(2))}.Location;
            Node_3 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(3))}.Location;
            
            Node_1_projected = PROJECTED_NODES(AT_NODES(DB_ELEMENTS{i}.Nodes(1)),:)';
            Node_2_projected = PROJECTED_NODES(AT_NODES(DB_ELEMENTS{i}.Nodes(2)),:)';
            Node_3_projected = PROJECTED_NODES(AT_NODES(DB_ELEMENTS{i}.Nodes(3)),:)';
            
            LOAD_FIELD(i,1+Vector_Field_Component) = Correction_Coefficient * Get_Tria_Element_Pressure(Node_1_projected, Node_2_projected, Node_3_projected, Pn, Span_length, Span_Pressure_Center, Root_chord_length, End_chord_length, Chord_Pressure_Center, A1, A2, A3, A4) * Get_Tria_Element_Area(Node_1_projected, Node_2_projected, Node_3_projected) / Get_Tria_Element_Area(Node_1, Node_2, Node_3);
            
        case 4
            
            Node_1 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(1))}.Location;
            Node_2 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(2))}.Location;
            Node_3 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(3))}.Location;
            Node_4 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(4))}.Location;
            
            Node_1_projected = PROJECTED_NODES(AT_NODES(DB_ELEMENTS{i}.Nodes(1)),:)';
            Node_2_projected = PROJECTED_NODES(AT_NODES(DB_ELEMENTS{i}.Nodes(2)),:)';
            Node_3_projected = PROJECTED_NODES(AT_NODES(DB_ELEMENTS{i}.Nodes(3)),:)';
            Node_4_projected = PROJECTED_NODES(AT_NODES(DB_ELEMENTS{i}.Nodes(4)),:)';
            
            LOAD_FIELD(i,1+Vector_Field_Component) = Correction_Coefficient * Get_Quad_Element_Pressure(Node_1_projected, Node_2_projected, Node_3_projected, Node_4_projected, Pn, Span_length, Span_Pressure_Center, Root_chord_length, End_chord_length, Chord_Pressure_Center, A1, A2, A3, A4) * Get_Quad_Element_Area(Node_1_projected, Node_2_projected, Node_3_projected, Node_4_projected) / Get_Quad_Element_Area(Node_1, Node_2, Node_3, Node_4);
            
        otherwise
            
            ErrCode = -1;
            fprintf (1, 'Недопустимое количество узлов [%i] у элемента с номером %i типа %s !\n', DB_ELEMENTS{i}.N_nodes, DB_ELEMENTS{i}.ID, DB_ELEMENTS{i}.type);
            break;
            
    end
        
end

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

fprintf (1, 'OK!\n');

Print_Delimiter();


% ======== Этап проверки нагружения ========

fprintf (1, 'ПРОВЕРКА НАГРУЖЕНИЯ\n');

% Создание системы координат для суммирования нагрузок

Set_Next_Coord_Number();

Summation_Coord = DB_Coord_Create_Rect_3P(Coord_Origin, Point_On_Axis_3, Point_On_Plane_13);

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

db_Next_Coord_Number = db_Next_Coord_Number + 1;

% Вычисление узловых сил

fprintf (1, '\nВычисление узловых сил...                       ');

NODE_FORCES = zeros(size(DB_NODES,1),1);

for i=1:size(DB_ELEMENTS,1)
    
    switch DB_ELEMENTS{i}.N_nodes
        
        case 3
            
            Node_1 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(1))}.Location;
            Node_2 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(2))}.Location;
            Node_3 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(3))}.Location;
            
            Node_1_projected = PROJECTED_NODES(AT_NODES(DB_ELEMENTS{i}.Nodes(1)),:)';
            Node_2_projected = PROJECTED_NODES(AT_NODES(DB_ELEMENTS{i}.Nodes(2)),:)';
            Node_3_projected = PROJECTED_NODES(AT_NODES(DB_ELEMENTS{i}.Nodes(3)),:)';
            
            Elm_Area = Get_Tria_Element_Area(Node_1, Node_2, Node_3);
            Elm_Pressure = LOAD_FIELD(i,1+Vector_Field_Component);
            Elm_Force = Elm_Pressure * Elm_Area;
            
            Mid_point = (Node_1_projected + Node_2_projected)/2;
            Centroid = Mid_point*2/3 + Node_3_projected/3;
            
            Node_Force_1 = Elm_Force * Get_Point_Line_Distance(Centroid, Node_2_projected, Node_3_projected) / Get_Point_Line_Distance(Node_1_projected, Node_2_projected, Node_3_projected);
            Node_Force_2 = Elm_Force * Get_Point_Line_Distance(Centroid, Node_3_projected, Node_1_projected) / Get_Point_Line_Distance(Node_2_projected, Node_3_projected, Node_1_projected);
            Node_Force_3 = Elm_Force - Node_Force_1 - Node_Force_2;
            
            NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(1))) = NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(1))) + Node_Force_1;
            NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(2))) = NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(2))) + Node_Force_2;
            NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(3))) = NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(3))) + Node_Force_3;
            
        case 4
            
            Node_1 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(1))}.Location;
            Node_2 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(2))}.Location;
            Node_3 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(3))}.Location;
            Node_4 = DB_NODES{AT_NODES(DB_ELEMENTS{i}.Nodes(4))}.Location;
            
            Elm_Area = Get_Quad_Element_Area(Node_1, Node_2, Node_3, Node_4);
            Elm_Pressure = LOAD_FIELD(i,1+Vector_Field_Component);
            Elm_Force = Elm_Pressure * Elm_Area;
            
            NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(1))) = NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(1))) + Elm_Force/4;
            NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(2))) = NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(2))) + Elm_Force/4;
            NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(3))) = NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(3))) + Elm_Force/4;
            NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(4))) = NODE_FORCES(AT_NODES(DB_ELEMENTS{i}.Nodes(4))) + Elm_Force/4;
            
        otherwise
            
            ErrCode = -1;
            fprintf (1, 'Недопустимое количество узлов [%i] у элемента с номером %i типа %s !\n', DB_ELEMENTS{i}.N_nodes, DB_ELEMENTS{i}.ID, DB_ELEMENTS{i}.type);
            break;
            
    end
    
end

Actual_Total_Force = sum(NODE_FORCES);

fprintf (1, 'OK!\n');

% Теоретический центр давления

Pressure_Center_XYZ_Projected = Get_Cartesian_Coord(A1, A2, A3, A4, Span_Pressure_Center, Chord_Pressure_Center);

% Фактический центр давления

fprintf (1, 'Вычисление фактического центра давления...      ');

Tmp_sum = zeros(3,1);

for i=1:size(NODE_FORCES,1)
    Tmp_sum = Tmp_sum + DB_NODES{i}.Location * NODE_FORCES(i);
end

Actual_Pressure_Center_XYZ = Tmp_sum / Actual_Total_Force;
Actual_Pressure_Center_XYZ_Projected = Project_Point_To_Plane(Actual_Pressure_Center_XYZ, Basic_Plane);
Actual_Pressure_Center_ksi = Get_Param_Coord(A1, A2, A3, A4, Actual_Pressure_Center_XYZ_Projected);
Actual_Pressure_Center_etta = Get_Param_Coord(A1, A3, A2, A4, Actual_Pressure_Center_XYZ_Projected);

fprintf (1, 'OK!\n');

% Фактические результирующие нагрузки

Actual_Total_Force_Vector = Actual_Total_Force * Load_Vector;
Actual_Total_Moment_Vector = cross(Actual_Pressure_Center_XYZ - Summation_Point, Actual_Total_Force_Vector);

FX = dot(Actual_Total_Force_Vector, DB_COORD{AT_COORD(Summation_Coord)}.e_x);
FY = dot(Actual_Total_Force_Vector, DB_COORD{AT_COORD(Summation_Coord)}.e_y);
FZ = dot(Actual_Total_Force_Vector, DB_COORD{AT_COORD(Summation_Coord)}.e_z);

MX = dot(Actual_Total_Moment_Vector, DB_COORD{AT_COORD(Summation_Coord)}.e_x);
MY = dot(Actual_Total_Moment_Vector, DB_COORD{AT_COORD(Summation_Coord)}.e_y);
MZ = dot(Actual_Total_Moment_Vector, DB_COORD{AT_COORD(Summation_Coord)}.e_z);

% Вывод статистики по нагрузкам

fprintf(1, '\nСуммарная сила:\n');
fprintf(1, '\n    Теоретическая: \t%.5f\n', Pn);
fprintf(1, '    Фактическая: \t%.5f\n', Actual_Total_Force);
fprintf(1, '    Отклонение: \t%.3f%%\n', (Actual_Total_Force/Pn - 1)*100);

fprintf(1, '\nЦентр давления:\n');
fprintf(1, '                                              Теор. знач.     Факт. знач.     Относ. откл.\n');
fprintf(1, '    Параметрическая координата по размаху:\t\t%.5f\t\t\t%.5f\t\t\t%.3f%%\n', Span_Pressure_Center, Actual_Pressure_Center_ksi, (Actual_Pressure_Center_ksi/Span_Pressure_Center - 1)*100);
fprintf(1, '    Параметрическая координата по хорде:\t\t%.5f\t\t\t%.5f\t\t\t%.3f%%\n', Chord_Pressure_Center, Actual_Pressure_Center_etta, (Actual_Pressure_Center_etta/Chord_Pressure_Center - 1)*100);

fprintf(1, '\n                                              Теор. знач.     Факт. знач.     Абсол. откл.\n');
fprintf(1, '    Проекция на базовую плоскость, X:\t\t\t%.5f\t\t%.5f\t\t%.5f\n', Pressure_Center_XYZ_Projected(1), Actual_Pressure_Center_XYZ_Projected(1), Pressure_Center_XYZ_Projected(1)-Actual_Pressure_Center_XYZ_Projected(1));
fprintf(1, '    Проекция на базовую плоскость, Y:\t\t\t%.5f\t\t%.5f\t\t%.5f\n', Pressure_Center_XYZ_Projected(2), Actual_Pressure_Center_XYZ_Projected(2), Pressure_Center_XYZ_Projected(2)-Actual_Pressure_Center_XYZ_Projected(2));
fprintf(1, '    Проекция на базовую плоскость, Z:\t\t\t%.5f\t\t%.5f\t\t%.5f\n', Pressure_Center_XYZ_Projected(3), Actual_Pressure_Center_XYZ_Projected(3), Pressure_Center_XYZ_Projected(3)-Actual_Pressure_Center_XYZ_Projected(3));

fprintf(1, '\n    Координаты в пространстве, [X,Y,Z]:    [%.5f, %.5f, %.5f]\n', Actual_Pressure_Center_XYZ(1), Actual_Pressure_Center_XYZ(2), Actual_Pressure_Center_XYZ(3));

fprintf(1, '\nФактические результирующие силовые факторы:\n');

fprintf(1, '\n    Полная сила: \t%.5f\n', Actual_Total_Force);
fprintf(1, '\n    Fx: \t%.5f\n', FX);
fprintf(1, '    Fy: \t%.5f\n', FY);
fprintf(1, '    Fz: \t%.5f\n', FZ);

fprintf(1, '\n    Полный момент: \t%.5f\n', norm(Actual_Total_Moment_Vector));
fprintf(1, '\n    Mx: \t%.5f\n', MX);
fprintf(1, '    My: \t%.5f\n', MY);
fprintf(1, '    Mz: \t%.5f\n', MZ);
    
Print_Delimiter();


% ======== Этап записи результатов ========

if (~strcmp(Output_File_Format, 'OFF'))

    fprintf (1, 'ЗАПИСЬ РЕЗУЛЬТАТОВ\n');

    switch Output_File_Format
        case 'CSV'
            Write_CSV(LOAD_FIELD, Work_Directory);
        case 'SES'
            Write_SES(LOAD_FIELD, Work_Directory, Field_Name);
        case 'BOTH'
            Write_CSV(LOAD_FIELD, Work_Directory);
            Write_SES(LOAD_FIELD, Work_Directory, Field_Name);
        otherwise
            ErrCode = -1;
            fprintf (1, 'Недопустимое значение параметра Output_File_Format [%s] !\n', Output_File_Format);
    end

    if (ErrCode < 0)
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
    end

    Print_Delimiter();

end

if (strcmp(Output_File_Format, 'OFF'))
    fprintf (1, 'Запись файла с полем нагрузок была отключена настройкой Output_File_Format = OFF\n');
end

fprintf (1, 'Работа программы завершена успешно!\n\n');

% Вывод времени выполнения программы
toc;

% Порядок переменных в функциях вычисления координат:

% ksi = Get_Param_Coord(A1, A2, A3, A4, X);
% etta = Get_Param_Coord(A1, A3, A2, A4, X);
% X = Get_Cartesian_Coord(A1, A2, A3, A4, ksi, etta)
